Vlasov model using kinetic phase point trajectories 
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A method of solution of the collisionless Vlasov equation, by following collisionless phase point 
trajectories in phase space, is presented. It is shown that by increasing the number of phase points, 
without enhancing the resolution of phase space grid, the accuracy of simulation will be improved. 
Besides, the phase points spacing introduces a smaller scale than grid spacing on which fine structures 
might be more conveniently handled. In order to perform simulation with a large population of phase 
points, an effective interpolation scheme is introduced that reduces the number of operations. It is 
shown that by randomizing initial position of the phase points along velocity axis, the recurrence 
effect does not happen. Finally, the standard problem of linear Landau damping will be examined. 
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There are two equally important numerical approaches 
to the kinetic problem. The first approach, particle-in- 
cell (PIC) methods, self-consistently models the plasma 
by a finite number of macroparticles on a fixed grid. Its 
advantage is replacing the Vlasov equation by the or- 
dinary differential equations of motions of macroparti- 
cles that makes PIC codes easily extendable to multi- 
dimensional applications. However, it is well-known that 
the numerical noise (proportional to 1/%/iV where N is 
the number of macroparticles) inherent to PIC simulation 
becomes, in some cases, too significant to allow a precise 
description of the distribution function (DF). Neverthe- 
less, the PIC approach can produce accurate results when 
a sufficiently large number of macroparticles involve in 
the simulation [1]. 

The second approach is direct solution of the Vlasov 
equation (for a collisionless plasma) that is noise-free in 
comparison to PIC simulation. The major problem of 
the Vlasov simulation has been the development of fine 
structures (filamentation) in velocity space; i.e., a prob- 
lem with no seemingly simple cure. Partial treatments 
such as increase in velocity resolution, have sharply lim- 
ited the ability to extend the above work to higher dimen- 
sions [2|44| and thus treat realistic problems. A large class 
of the Vlasov simulation models is based on discretizing 
the Vlasov equation (mainly by a splitting scheme) on a 
phase space fixed grid. In the splitting method, the new 
/ was obtained as an algebraic expression in terms of 
the old / by a suitable interpolation method ||. There 
are a couple of problems with the splitting scheme: (i) 
it was not rigorously shown under which circumstances 
the coupled equations have solutions approximately con- 
sistent with the Vlasov equation; and (ii) following char- 
acteristics along the phase space coordinates departs one 
from the characteristics on which / truly remains invari- 
ant [3(| . Besides, for them the partial treatment of the 
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filamentation, i.e. increasing the velocity resolution, is 
very expensive because all of the simulation operations 
is performed on the grid. Moreover, they are suffering 
from the recurrence effect [6j]. 

Integration of the Vlasov equation along the collision- 
less phase point trajectories has been the most promising 
of these methods Q. This method is based on following 
the characteristics along which / is constant in the col- 
lisionless case. Therefore, characteristic equations of the 
Vlasov equation are solved. That is one of the advan- 
tages of the scheme and makes it possible to use all the 
PIC simulation experiences. A complete description of 
the method of characteristics was presented in Ref. Q . 
There, besides the grid points there are equal number of 
phase points over which the DF, f p , is initially defined. 
Interpolation is performed between the phase points and 
the fixed background grid to obtain the DF, f g , on the 
grid. The main advantage of introducing f p is that it 
remains unchanged when the phase points follow their 
characteristics (contrary to the semi-Lagrangian method 
in which the DF is altered by the interpolation). 

In present paper, we first improve the accuracy of the 
simulation by increasing the number of phase points in 
each grid cell, without enhancing the resolution of the 
phase space grid. Better accuracy is the result of sam- 
pling DF with higher resolution. Larger population of 
phase points, in comparison to grid points, introduces a 
smaller scale than grid spacing on which fine structures 
might be more conveniently handled. However, increas- 
ing the number of phase points necessitates an effective 
interpolation scheme (IS) that reduces the number of op- 
erations while keeping the accuracy [e.g. in comparison 
to the bilinear interpolation scheme (BIS)]. Accordingly, 
a new IS is introduced. In such a way that the DF of 
each grid point is obtained by averaging the DF of phase 
points located in four cells around a grid point. More- 
over, it is shown that by randomizing the initial position 
of the phase points along velocity axis, the recurrence 
effect does not happen and the reason is given in detail. 

Our mathematical model is the one-dimensional 
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FIG. 1: A typical part of the phase space grid, (a) Regu- 
lar arrangement (b) Random arrangement along the velocity 
direction. 
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FIG. 2: The simulation result of the free streaming part of 
the Vlasov equation with BIS. (a) Regular arrangement (b) 
Random arrangement along the velocity direction. 



Vlasov-Poisson system, 

d t f + vd x f-E{x,t)d v f = 0, 

/+oo 
fdv, (1) 
-oo 

where / is the electron DF and E is the electric field. In 
Eq. (1), and in the rest of the article, time is normalized 
to the inverse electron plasma frequency ui~ e , space is 
normalized to the Debye length Ad, and velocity is nor- 
malized to the electron thermal speed vt £ = ^DW pe . Ions 
are taken to be motionless, and their only role is to pro- 
vide a uniform, neutralizing background. Furthermore, 
periodic boundary conditions are assumed in x. 

As it was mentioned, the present solution of Vlasov 
equation is based on following the phase points trajec- 
tories along which phase-space DF is constant. In order 
to obtain the phase point trajectory one has to solve the 
characteristics of the Vlasov equation, 

dxp/dt = v p , 

dvp/dt = ~E p , (2) 

where subscript "p" stands for "phase point" . 

Let us first begin with the free streaming part 
(the advection term) of the Vlasov equation, dtf + 
vd x f — 0, through the following example. The so- 
lution of the advection part at a time t is given 
as a function of the initial condition by the relation 
f(x,v,t) = f(x — vt,v,0). If we consider an initial 
Maxwellian distribution perturbed by a small pertur- 
bation, f(x, v, 0) = l/v / 2~7rexp(— v 2 /2) [1 + ecos(fcx)], 
then the charge density (J_°° fdv — 1) will be given by 

p(x,t) = ex-p(—k 2 t 2 )ecos(kx) [6(. The analytical solu- 
tion is decaying exponentially in time. For the numerical 
solution, the charge density is calculated at every spatial 
grid point by summation over all grid points in velocity 
space. On the Eulerian grid due to equal spacing along 
the velocity axis Av, the initial condition can be recon- 
structed at recurrence time, Tr = 2Tr/(kAv) @. 

Now, let us examine the above analytical solution by 
the method of characteristics. For the free streaming, we 
just need to solve dx/dt = v. Since in free streaming the 



velocity of the phase points is constant, the characteristic 
equation can be exactly solved, that is = £™ + v p t n , 
where the superscript denotes t = nAt. We first con- 
sider a fixed grid with regular phase points arrangement 
(Fig. la). Second, according to each x p and v p , its f p is 
allocated. Next, x p is advanced one time step while v p 
remains constant. Then, interpolation is performed be- 
tween the phase points and the fixed grid in phase space 
by BIS [Ho} to obtain f g . Finally, the charge density is 
obtained by summation (here, Trapezoidal rule) over all 
grid points in velocity space. 

To compare the simulation result, we choose the pa- 
rameters similar to Ref. That is, e = 0.1, with grid 
points N x x N v — 16 x 32. The length in space is L = An 
and in velocity space we use — 5 < v < 5. The recurrence 
time of the Eulerian codes, for this case, is Tr = 38.95. 
We put hundred phase points (10 x 10) in each grid cell. 
The first result of the model was surprising. The recur- 
rence took place at 442.3362 instead of 38.95 (Fig. 2a). 
We realized that contrary to the Eulerian codes, it is not 
the velocity grid that specifics the recurrence time. 

Thus, the main question is "How should the recurrence 
time be calculated?" To answer this question we have to 
note that it is the evolution of the phase points arrange- 
ment that changes the interpolation weighting and there- 
fore f g (also p). Thus, as time goes on, the recurrence 
will take place if there is a possibility to reconstruct the 
initial phase points arrangement (arrangement at t = 0). 
According to Fig. la, the velocity spacing in regular 
arrangement is dpv. That means, the smallest velocity 
for the moving phase point is dpv (in the positive direc- 
tion) and the other phase points velocities are integer 
multiples of dpv. Since the boundary condition is pe- 
riodic, if those phase points that their velocity is dpv , 
move a distance L within a time interval Tr, the other 
phase points will move an integer multiples of L within 
the same time interval. As a result, all the phase points 
return to their initial positions and the recurrence oc- 
curs. Accordingly, in our model, the recurrence time is 
Tr = L/dpv. In another word, there is a much smaller 
scale dpv, in comparison to the grid spacing, along the 
velocity axis that becomes an important factor in the dy- 
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FIG. 3: The simulation result of the free streaming part of 
the Vlasov equation with AIS. (a) Regular arrangement (b) 
Random arrangement along the velocity direction. 



namics. The influence of dpv on the recurrence time is an 
evidence regarding our claim that fine structures might 
be more conveniently handled. For the above example, 
dpv = 0.0284 that leads to correct recurrence time (Fig. 
2a). It is obvious from the Fig. 2a that there are several 
other recurrences with smaller amplitude. According to 
our analysis, these smaller amplitude recurrences are due 
to sub-arrangements of the phase points, that is, a small 
group of phase points (not all) has return to their original 
arrangement. 

The latter analysis of the recurrence time is based on 
two facts. First, the boundary condition is periodic. Sec- 
ond, the velocity of each phase point is an integer multi- 
ples of the others. Therefore, by randomizing the phase 
points velocities (Fig. lb), we can prevent the occurrence 
of recurrence (Fig. 2b). To do that, we used a random 
generator to modify the velocity of each phase point in 
the range [— dpv/2, dpv/2] (Fig. lb). Figure 2 depicts 
the result after randomizing the velocities of the phase 
points. The recurrence that is supposed to take place at 
442.3362 does not happen. That means our suggestion 
for the recurrence mechanism works. 

As it was mentioned, increasing the number of phase 
points necessitates an effective IS that reduces the num- 
ber of operations while keeping the accuracy (e.g. in 
comparison to BIS). For this purpose, we develop a new 
IS that is almost as accurate as BIS with the difference 
that it does not use the weighting mechanism. 

Let us call the new IS as average interpolation scheme 
(AIS). In AIS like BIS, all the phase points should be 
swept one by one to find those that are located in four 
cells around a specific grid point (with the important 
difference that in AIS, finding the host cell of the phase 
point is the only step, please see below). Now, we denote 
the distance between the ith phase point (in those four 
cells) and the grid point along x axis by Axi and along 
v axis by A«j. Then, by using the Taylor expansion of 



FIG. 4: (a) Comparison of BIS and AIS. (b) Improvement of 
the accuracy by increasing the number of phase points. 



all J phase points, located in four cell, we obtain 
1 J 

fg( x g> v g) = f P t(x g + Ax u v g + A Vl ) 

J i=i 

-(9xfg)j E A ^ - ^h)! E A ^ + °( A2 )- 



i=l 



For an optimum total number of phase points, their ini- 
tial density in phase space is almost uniform and there- 
fore, J2i=i — an d J2i=i — 0- Moreover, the 
factor of 1/ J makes the approximation better. According 
to Liouville's theorem the density of system points in the 
vicinity of a given system point traveling through phase- 
space is constant with time. Therefore, the uniformity of 
phase points density is almost a constant of motion (to 
the extent of the truncation error). That means, 



fg ( x g ) v g ) 



1 J 



0(A 2 ). 



fpi (, x g ~t~ Axi 



-Avi) around (x g , v g ) and summing over 



Note, the phase points positions do not explicitly inter- 
fere in AIS (as through weighting comes to play in BIS). 
This is an essential feature of AIS that first makes AIS 
easily extendable to higher dimensions and second re- 
duces the number of operations. 

Figure 3 demonstrates the simulation of free stream- 
ing that is performed using AIS. Fig. 3a shows the result 
when the initial phase points arrangement is regular. At 
it is expected a recurrence happen at t = 442.3362. How- 
ever, after randomizing the initial phase point arrange- 
ments, the recurrence does not take place. 

Now, we first compare the accuracy of BIS and AIS and 
then show that the accuracy of the method will improve 
by increasing the number of phase points. In order to 
compare BIS and AIS, we perform two simulations with 
these ISs. Both of simulations are initially fed by ran- 
dom phase points arrangement with 2500 phase points 
in each cell. Fig. 4(a) exhibits the result. Although, BIS 
is slightly more accurate (maximum 1 x 10 -4 ) but AIS 
interpolates much faster than BIS. It is clear that after 
initial stage the accuracy of two interpolations becomes 
almost similar. Note that in this letter we use the sim- 
plest scheme with the parameters that are not necessarily 
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FIG. 5: The linear Landau damping (a) The exponential de- 
cay of the amplitude of the electric field with the damping 
rate of 7 = 0.153 (b) The relative error in the total energy in 
percent. 

the most appropriate ones (Ax = 0.83 and Av = 0.3125. 
Recall that the accuracy of the interpolation schemes is 
0(A 2 ). Therefore, we have not invested on the accu- 
racy and the parameters have chosen to make our results 
comparable with the results of Ref. Q. Next, in order 
to show the improvement of accuracy of the method, we 
redo the simulation of free-streaming by AIS for two dif- 
ferent cases when 100 and 2500 phase points are put in 
each cell. Figure 4(b) demonstrates the results. It is ob- 
vious that when we put 2500 phase points in each cell, 
the result is one order of magnitude more accurate than 
the case when the DF is sampled by 100 phase points 
within a cell. 

We now examine classical numerical test of the linear 
Landau damping. In this case, we start with /(x, v, 0) = 
l/v / 27rexp(— v 2 /2) [1 + ecos(fcx)]. The initial arrange- 
ment of the phase points is random along w-axis and 
regular along x-axis. By AIS, f g is calculated. In- 
tegrating f g over the velocity space and solving Pois- 
son's equation leads to the electric field on the grid, E g . 
Poisson's equation is solved by the Fast Fourier Trans- 
form. Then, by a Lagrange polynomial IS Q, E g is 



interpolated to the position of phase points to obtain 
E p . Having E p , Eqs. (2) are solved by the Leapfrog- 
Trapezoidal scheme. The parameters are e = 0.01, with 
grid points N x x N v = 256 x 256 and in each cell we put 
4 phase points. The periodic length is L — 4ir, k = 0.5, 
—5 < v < 5 , and dt = 0.1. In Fig. 5a, the basic mode 
of the electric field is plotted against time. It shows the 
exponential decay of the amplitude of the electric field ac- 
cording to Landau's theory. The damping rate, the slope 
of straight line, obtained by this method is 7 = 0.153 
which agrees very well with values predicted by the the- 
ory Figure 5b exhibits the relative error in the total 
energy in percent, i.e. (total energy"-total energy ) /total 
energy xlOO, recalling that the superscript "ra" denotes 
the quantities at t — nAt). As it is seen, the scheme is 
capable of keeping the energy conservation. 

In summary, we improved the method of characteris- 
tics [H by increasing the number of phase points in phase 
space. Naturally, due to phase points dynamics the de- 
velopment of steep gradients can be managed without 
enhancing the velocity resolution. In order to reduce the 
number of operation AIS is introduced that is easily ex- 
tendable to higher dimension. By AIS and randomizing 
the phase point arrangement, we could prevent the oc- 
currence of the recurrence effect. The scheme is very 
similar to PIC simulation with the advantages that it is 
noise-free with a simpler IS (i.e. weighting is omitted) 
and that through DF it is much easier to calculate the 
thermodynamic quantities. Using todays supercomput- 
ers, this method appears to be a good alternative to the 
PIC methods for dealing with strongly nonlinear prob- 
lems in phase space when little noise and good precision 
is needed. 
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